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Abstract. The syntactic nature and compositionality characteristic of stochastic process 
algebras make models to be easily understood by human beings, but not convenient for 
machines as well as people to directly carry out mathematical analysis and stochastic sim- 
ulation. This paper presents a numerical representation schema for the stochastic process 
algebra PEPA, which can provide a platform to directly and conveniently employ a variety 
of computational approaches to both qualitatively and quantitatively analyse the models. 
Moreover, these approaches developed on the basis of the schema are demonstrated and 

(N discussed. In particular, algorithms for automatically deriving the schema from a general 

O PEPA model and simulating the model based on the derived schema to derive performance 

g measures are presented. 
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T — I 

r— I 1 Introduction 

o 

Stochastic process algebras, such as PEPA [l], TIPP [2], EMPA [s], are powerful modelling for- 
ce malisms for concurrent systems which have enjoyed considerable success over the last decade. 

A stochastic process algebra model is constructed to approximately and abstractly represent a 
^ system whilst hiding its implementation details. Based on the model, performance properties 

^ of the dynamic behaviour of the system can be assessed, through some techniques and com- 

putational methods. This process is referred to as the performance modelling of the system, 
which mainly involves three levels: model construction, technical computation and performance 
derivation. In order to derive performance measures from large scale stochastic process algebra 
models, many mathematical tools and approaches have been proposed to study the models. 
T— I For instance, a fluid approximation method has been proposed in f4] to avoid the state space 

explosion problem encountered in the analysis of large scale PEPA models. 
• • However, the syntactic nature of stochastic process algebras makes models easily understood 

by human beings, but not convenient for machines/computers (as well as for human beings) to 
^ directly employ these tools and approaches. In addition, the compositionality of the formalisms 

allow a model to be locally defined, but the analysis of the model or the underlying continuous- 
time Markov chain (CTMC) carries out in the global manner since it is the whole system rather 
than a part of it to be usually interested and considered. The syntactical and compositional 
qualities of the stochastic process algebras, which are advantages in model construction, turn 
to be disadvantages in model analysis. 



* School of Information Engineering, Yangzhou University, Yangzhou 225009, P. R. China. Email: jied- 
ing@yzu.edu.cn 

^Laboratory for Foundations of Computer Science, School of Informatics, The University of Edinburgh, Ed- 
inburgh, UK. E-mail: jane.hillston@ed.ac.uk 



1 



1.1 Paper contributions 



In order to overcome the obstacles in the direct and convenient apphcation of the mathematical 
tools, we propose a new numerical representation schema for the formulism PEPA in this paper. 
In this schema, labelled activities are defined to cope with the difference between actions in 
PEPA and transitions in the underlying CTMC, so that the correspondence between them is 
one-to-one. Activity matrices based on the labelled activities are defined to capture structural 
information about PEPA models. Moreover, transition rate functions are proposed to capture 
the timing information. These concepts numerically describe and represent a PEPA model, 
and provide a platform for conveniently and easily simulating the underlying CTMC, deriving 
the fluid approximation, as well as leading to an underlying Place/Transition (P/T) structure. 
These definitions are consistent with the original semantics of PEPA, and a PEPA model can 
be recovered from its numerical representation. An algorithm for automatically deriving the 
schema from any given PEPA model has been provided. Some characteristics of this numerical 
representation are revealed. For example, using numerical vector forms the exponential increase 
of the size of the state space with the number of components can be reduced to at most a 
polynomial increase. 




Figure 1: Three levels of PEPA modelling 

The benefits of the schema embodies the three aspects of performance modelling, which are 
illustrated by Figure [T} At the first level, the proposed new schema numerically describes any 
given PEPA model and provides a platform to directly employ a variety of approaches to analyse 
the model. These approaches are shown at the second level. At this level, a fluid approximation 
method for the quantitative analysis of PEPA is established, as well as investigated, mainly with 
respect to its convergence and the consistency between this method and the underlying CTMC. 
In addition, a P/T structure-based approach is revealed, which can be utilised to qualitatively 
analyse the model. At the third level, both qualitative and quantitative performance measures 
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can be derived from the model through those approaches. A stochastic simulation algorithm 
for the aggregated CTMC, which is based on the numerical representation schema, is proposed 
to obtain general performance metrics in this paper. As for the other two approaches, related 
investigation and analysis were given in [5] and [g] respectively, which will be briefly introduced 
in the next subsection. 

1.2 Related work 

Our work is motivated and stimulated by the pioneering work on the numerical vector form 
and activity matrix in [1], which was dedicated to the fluid approximation for PEPA. The P/T 
structure underlying each PEPA model, as stated in Theorem [T] in this paper, reveals tight con- 
nections between stochastic process algebras and stochastic Petri nets. Based on this structure 
and the theories developed for Petri nets, several powerful techniques for structural analysis of 
PEPA were presented in 1 5 ] , including a structure-based deadlock-checking method which avoids 
the state space explosion problem. In [T], a new operational semantics was proposed to give 
a compact symbolic representation of PEPA models. This semantics extends the application 
scope of the fluid approximation of PEPA by incorporating all the operators of the language 
and removing earlier assumptions on the syntactical structure of the models amenable to this 
analysis. Moreover, the paper |6| shows how to derive the performance metrics such as action 
throughput and capacity utilisaition from the fluid approximation of a PEPA model. 

1.3 Paper organisation 

The remainder of this paper is structured as follows: Section 2 gives a brief introduction to the 
PEPA formulism; In Section 3, 4 and 5, we respectively present the three combinators of the 
numerical schema, i.e. the numerical vector form, labelled activity and activity matrix, as well 
as the transition rate function. Computational approaches for performance derivation that are 
developed on the basis of the schema are demonstrated in Section 6. We flnally conclude the 
paper in Section 7. 

2 Introduction to PEPA 

PEPA (Performance Evaluation Process Algebra) [l], developed by Hillston in the 1990s, is a 
high-level model speciflcation language for low-level stochastic models, and describes a system 
as an interaction of the components which engage in activities. In contrast to classical process 
algebras, activities are assumed to have a duration which is a random variable governed by an 
exponential distribution. Thus each activity in PEPA is a pair (a, r) where a is the action 
type and r is the activity rate. The language has a small number of combinators, for which we 
provide a brief introduction below; the structured operational semantics can be found in 
The grammar is as follows: 

S ::= {a,r).S\S + S\Cs 
P ::= PtXiP \ P/L \ C 

where S denotes a sequential component and P denotes a model component which executes 
in parallel. C stands for a constant which denotes either a sequential component or a model 
component as introduced by a deflnition. Cg stands for constants which denote sequential 
components. The effect of this syntactic separation between these types of constants is to 
constrain legal PEPA components to be cooperations of sequential processes. 
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Prefix: The prefix component {a,r).S has a designated first activity (a,r), which has 
action type a and a duration which satisfies exponential distribution with parameter r, and 
subsequently behaves as S. 

Choice: The component S + T represents a system which may behave either as S or as 
T. The activities of both S and T are enabled. Since each has an associated rate there is a 
race condition between them and the first to complete is selected. This gives rise to an implicit 
probabilistic choice between actions dependent of the relative values of their rates. 

Hiding: Hiding provides type abstraction, but note that the duration of the activity is 
unaffected. In P/L all activities whose action types are in L appear as the "private" type r. 

Cooperation:P IXl Q denotes cooperation between P and Q over action types in the coop- 
eration set L. The cooperands are forced to synchronise on action types in L while they can 
proceed independently and concurrently with other enabled activities {individual activities). 
The rate of the synchronised or shared activity is determined by the slower cooperation (see [T] 
for details). We write P \\ Q as an abbreviation for P 1X1 Q when L = ^ and P[N] is used to 
represent N copies of P in a parallel, i.e. P[3] = P || P || P. 

Constant: The meaning of a constant is given by a defining equation such as A = P. This 
allows infinite behaviour over finite states to be defined via mutually recursive definitions. 

On the basis of the operational semantic rules (please refer to [l] for details) , a PEPA model 
may be regarded as a labelled multi-transition system 



where C is the set of components, Act is the set of activities and the multi-relation is given 
by the rules. If a component P behaves as Q after it completes activity (a,r), then denote the 

transition as P ^——^ Q. 

The memoryless property of the exponential distribution, which is satisfied by the durations 
of all activities, means that the stochastic process underlying the labelled transition system 
has the Markov property. Hence the underlying stochastic process is a CTMC. Note that in 
this representation the states of the system are the syntactic terms derived by the operational 
semantics. Once constructed the CTMC can be used to find steady-state or transient probability 
distributions from which quantitative performance can be derived. 

3 Numerical Vector Form 

The usual state representation in PEPA models is in terms of the syntactic forms of the model 
expression. When a large number of repeated components are involved in a system, the state 
space of the CTMC underling the model can be large. This is mainly because each copy of the 
same type of component is considered to be distinct, resulting in distinct Markovian states. The 
multiple states within the model that exhibit the same behaviour can be aggregated to reduce 
the size of the state space as shown by Gilmore et al. (s] using the technique based on a vector 
form. The CTMC is therefore constructed in terms of equivalence classes of syntactic terms. 
"At the heart of this technique is the use of a canonical state vector to capture the syntactic 
form of a model expression", as indicated in (4j. Rather than the canonical representation 
style, an alternative numerical vector form was proposed by Hillston in [i] for capturing the 
state information of models with repeated components. In the numerical vector form, there 
is one entry for each local derivative of each type of component in the model. The entries in 
the vector are the number of components currently exhibiting this local derivative, no longer 
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syntactic terms representing the local derivative of the sequential component. Following [Ij, 
hereafter the term local derivative refers to the local state of a single sequential component. 

Definition 1. (Numerical Vector Form [41)- For an arbitrary PEPA model M. with n 
component types Ci,i = 1,2, • • • ,n, each with di distinct local derivatives, the numerical vector 
form of M, m(A^), is a vector with d = X^iLi di entries. The entry m[Ci ] records how many 
instances of the jth local derivative of component type Ci are exhibited in the current state. 

By adopting this model-aggregation technique, the number of the states of the system can be 
reduced to only increase (at most) polynomially with the number of instances of the components. 
According to Definition [l| m{Ci-) > for each d-. At any time, each sequential component 
stays in one and only one local derivative. So the sum of m(Ci), i.e. Y^'jLi ^[Cij], specifies the 
population of Ci in the system. Notice that m(Cj) satisfies 



EU^[Ci,] = Mi. 



(3.1) 



Then according to the well-known combinatorial formula (Theorem 3.5.1 in [o]), there are 

( ^ ^\ ^1 solutions, i.e. ( ^ ^1 states in terms of C,- in the system. But the 
\ di-1 J ' \ di-1 J 

possible synchronisations in the PEPA model have not been taken into account in the restrictions 
(3.1) and thus the current restrictions may allow extra freedom for the solutions, so the given 

number ^ ^ ^ is an upper bound of the exact number of the states in terms of Ci. 

Notice that 

+ -1 \ (Mi + di-1)! 



di-1 J (di- l)!Mi! 

< (M, + di- . 

Therefore, it is easy to verify the following 

Proposition 1. Consider a system comprising n types of component, namely Ci,C2,-- - ,Cn, 
with Mi copies of the component of type Ci in the system, where Ci has di local derivatives, for 
i = 1,2, ■ ■ ■ ,n. Then the size of the state space of the system is at most 



i=i ^ * ^ i=i 



The upper bound given in Proposition [T] guarantees that the size of the state space increases 
at most polynomially with the number of instances of the components. Consider the following 
PEPA model. 



Model 1. 



User I ={taski,a).U ser2 
User2 ={task2,b).U seri 
Severi ={taski,a).Sever2 
Sever2 ={reset, d). Severi 
UserilM] ^^^^^^ Severi[N]. 
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(2,0,0,2)' 




(1,1,1,1)^ 



0,2,0,2) 



(0,2,1,1) 



(0,2,2,0)'' 



Figure 2: Transition between states (a revised version of the one in 

According to the semantics of PEPA originahy defined in [I], the size of the state space of 
the CTMC underlying Model [l] is 2 ^ , which increases exponentially with the numbers of the 
users and severs in the system. According to Definition [T| the system vector m has four entries 
representing the instances of components in the total four local derivatives, that is 



m = {in.[Useri],in[User2],in.[Severi],ni[Sever2]) ■ 

Let M = N = 2, then the system equation of Model [l] determines the starting state (2, 0, 2, 0). 
By enabling activities or transitions, all reachable system states can be manifested, see Figure[2] 
The size of the state space is nine. The upper bound of the size given by Proposition [T| 

^ ^ ^ ^)(^^2~^1 ^^"^^ ("^^ + 1) X (N + 1), is nine, coinciding with the size of the 
state space. The bound given in Proposition [l] is sharp and can be hit in some situations. 



4 Labelled Activity and Activity Matrix 

In the PEPA language, the transition is embodied in the syntactical definition of activities, in 
the context of sequential components. Since the consideration is in terms of the whole system 
rather than sequential components, the transition between these system states should be defined 
and represented. This section presents a numerical representation for the transitions between 
system states and demonstrates how to derive this representation from a general PEPA model. 

If a system vector changes into another vector after firing an activity, then the difference 
between these two vectors manifests the transition corresponding to this activity. Obviously, 
the difference is in numerical forms since all states are numerical vectors. Consider Model [l] 
and its transition diagram in Figure [2} Each activity in the model corresponds to a vector, 
called the transition vector. For example, taski corresponds to the transition vector = 
(— 1, 1, —1, 1)-^. That is, the derived state vector by firing taski from a state, can be represented 
by the sum of l^"-^'^^ and the state enabling taski. For instance, (2, 0, 2, 0)^ + /*"*'^i = (1, 1, 1, 1)-^ 
illustrates that si = (2, 0, 2, 0)^ transitions into S2 = (1, 1, 1, 1)^ after enabling taski. 

Similarly, task2 corresponds to l^""^'^^ = (1,-1,0,0)"^ while reset corresponds to T'^'*^* = 
(0,0, —1, 1). The three transition vectors form a matrix, called the activity matrix, see Table [ij 
Each activity in the model is represented by a transition vector — a column of the activity 
matrix, and each column expresses an activity. So the activity matrix is essentially indicating 
both an injection and a surjection from syntactic to numerical representation of the transition 
between system states. The concept of the activity matrix for PEPA was first proposed by Hill- 
ston in [4l[l0|. However, the original definition cannot fully reflect the representation mapping 
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Table 1: Transition vectors form an activity matrix 





jtaski 


jtask2 


jreset 




Useri 


-1 


1 







User2 


1 


-1 







Severi 


-1 





1 




Sever2 


1 





-1 





considered here. This is due to the fact of that the original definition is local-derivative-centric 
rather than transition centric. This results in some limitations for more general applications. 
For example, for some PEPA models (e.g. Model [2] in the following context), some columns of 
the originally defined matrix cannot be taken as transition vectors so that this definition cannot 
fully reflect the PEPA semantics in some circumstances. In the following, a modified definition 
of the activity matrix is given. The new definition is activity- or transition-centric, which brings 
the benefit that each transition is represented by a column of the matrix and vice versa. 

Model 2. 

P/=VO-^2 + (a,r::).P3 

P2={P,rp).P, + {P,r'p).P^ 

P3=(7,r,).Pi 

Qi ={a,ra_).Q2 

Q2 ={-f,r'^).Qi 

P,[A]^MB]. 




Figure 3: Transition diagram of Model [2] 

Let us consider a PEPA model, i.e. the following Model [2| in which there are multiple choices 
after firing some activities. In this model, firing a in the component P may lead to two possible 
local derivatives: P2 and P3, while firing f5 may lead to Pi and P3. In addition, firing 7 may lead 
to Pi, Qi. However, only one derivative can be chosen after each firing of an activity, according 
to the semantics of PEPA. But the original definition of activity matrix cannot clearly reflect 
this point. See the activity matrix of Model[2]given in Table[2j Moreover, the individual activity 
7 in this table, which can be enabled by both P3 and Q2, may be confused as a shared activity. 

In order to better reflect the semantics of PEPA, we modify the deflnition of the activity 
matrix in this way: if there are m possible outputs, namely {Pi, P2, • • • , Rm}, after firing either 
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Table 2: Originally defined activity matrix of Model [2] 





a 


/3 


7 


Pi 


-1 


1 


1 




1 


-1 





P3 


1 


1 


-1 


Qi 


-1 





1 


Q2 


1 





-1 



an individual or a shared activity /, then I is "split" into m labelled Is: I'^^^V"'^^ ■ ■ ■ Here 
{wi}^^ are m distinct labels, corresponding to {Ri}^^ respectively. Each l"^' can only lead to 
a unique output Ri. Here there are no new activities created, since we just attach labels to the 
activity to distinguish the outputs of firing this activity. The modified activity matrix clearly 
reflects that only one, not two or more, result can be obtained from firing /. And thus, each 
can represent a transition vector. 

For example, see the modified activity matrix of Model [2] in Table [3} In this activity matrix, 
the individual activity 7 has different "names" for different component types, so that it is not 
confused with a shared activity. Another activity /3, is labelled as /3^2-s>Pi ^j^^ pP2^P3^ 
respectively reflect the corresponding two choices. In this table, the activity a is also split and 
attached with labels. 



Table 3: Modified activity matrix of Model [2] 





t 
r 

a 


CM 

t 

T 
a 
















^P3-^Pl 




Pi 


-1 


-1 


1 





1 





P2 


1 





-1 


-1 








P3 





1 





1 


-1 





Qi 


-1 


-1 











1 


Q2 


1 


1 











-1 



Before giving the modified definition of activity matrix for any general PEPA model, the 
pre and post sets for an activity are first defined. For convenience, throughout this paper any 

{l,r) (l^rf^^) 

transition U -—4 V defined in the PEPA models may be rewritten as U — > V, or just 
U — ^ V if the rate is not being considered, where U and V are two local derivatives. 

Definition 2. (Pre and post local derivative) 

1. If a local derivative U can enable an activity I, that is U — ^ •, then U is called a pre local 
derivative of I. The set of all pre local derivatives of I is denoted by pre(/), called the pre 
set of I. 

2. If V is a local derivative obtained by firing an activity I, i.e. -—^^V, then V is called a 
post local derivative of I. The set of all post local derivatives is denoted by post(/), called 
the post set of I. 
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3. The set of all the local derivatives derived from U by firing I, i.e. 



post(C/,/) = {V \ U ^V}, 

is called the post set of I from U . 

Obviously, if / has only one pre local derivative, i.e. #pre(/) = 1, then / is an individual 
activity, like /3 in Model [2j whereafter the notation is defined as the cardinality of the set 
A, i.e. the number of elements of A. But I being individual does not imply #pre(/) = 1, see 7 
for instance. If I is shared, then #pre(Z) > 1, for example, see #pre(a) = Qi} = 2. For a 

shared activity I with pre(/) = k, there are k local derivatives that can enable this activity, each 
of them belonging to a distinct component type. The obtained local derivatives are in the set 
post(pre(/)[i], Z), where pre(/)[i] is the i-th pre local derivative of /. But only one of them can 
be chosen after I is fired from pre(/)[i]. Since for the component type, namely i or Cj, there are 
#post(pre(/)[i], outputs, so the total number of the distinct transitions for the whole system 
is nf=i #P°s^(P^6(0 0- That is, there are 11^=1 #Post(pre(/)[z], /) possible results but only 
one of them can be chosen by the system after the shared activity / is fired. In other words, to 
distinguish these possible transitions, we need Y[i=i 7^Post(pre(/)[i], I) different labels. Here are 
the readily accessible labels: 

(pre(0[l] ^ Fi,pre(Z)[2] ^ Fs,--- ,pre(/)[A:] ^ 

where Vi S post(pre(^)[2], Obviously, for each vector (Vi,V2,-- - ,Vfc) in post(pre(/)[l], Z) x 
post(pre(/)[2],/)x---xpost(pre(/)[A;],/),the labelled activity /(P''«=WW^^i'P'^^(')[2Hy2,- ,prc(0[fe]^v,) 
represents a distinct transition. For example, a in Model [2] can be labelled as q;(^i~^^2'Qi~''Q2) 
and a(^i^^3,Qi^Q2). 

For an individual activity I, things are rather simple and easy: for U € pre(Z), / can be 
labelled as /^^post(f/,0[i]^ iU^post(um ^ ^c/^post(r/,z)[%]^ ^^^^^ ^ #post(;7,/). Varying U G 

pre(Z), there are Yluepre{i) ?^POst(C/, /) labels needed to distinguish the possible transitions. See 
pP2^Pi ^ l]P2^P3 ^ryPs^Pi ^jQ2^Qi Model ^ for instance. Now we give the formal definition. 



Definition 3. (Labelled Activity). 

1. For any individual activity I, for each U G pre(Z),y G post(f/, /), label I as 

2. For a shared activity I, for each 

(Vi,V2,-- - ,Vfc) G post(pre(0[l],0 x post(pre(0[2], x ••• x post(pre(0[fc], 0, 

label I as V" , where 

u; = (pre(0[l] ^Fi,pre(0[2] ^Fa,--- ,pre(OW ^Ffc). 

Each l^^^ or l^ is called a labelled activity. The set of all labelled activities is denoted by 
-^labei- For the above labelled activities l^^^ and l"^ , their respective pre and post sets are 
defined as 

pre^^-^) = {t/}, post(/^-^) = {y}, 
pre(r) = pre(0, post(r) = {Vi,V2,--- , Vk}. 

According to Definition [sj each l^^^ or l"^ can only lead to a unique output. No new 
activities are created, since labels are only attached to the activity to distinguish the results 
after this activity is fired. 

The impact of labelled activities on local derivatives can be recorded in a matrix form, as 
defined below. 
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Definition 4. (Activity Matrix). For a model with Nj^^^^^^ labelled activities and Njy distinct 
local derivatives, the activity matrix C is an Nx> x NyXi^^^i matrix, and the entries are defined 
as follows 

{+1 if Ui £ post{lj) 
-1 ^fU^£pre{lJ) 
otherwise 

where Ij is a labelled activity. The pre activity matrix C^''^ and post activity matrix CP°®* are 
defined as 

The modified activity matrix captmes ah the structural information, including the informa- 
tion about choices and synchronisations, of a given PEPA model. From each row of the matrix, 
which corresponds to each local derivative, we can know which activities this local derivative 
can enable and after which activities are fired this local derivative can be derived. From the 
perspective of the columns, the number of "— l"s in a column tells whether the corresponding 
activity is synchronised or not. Only one "—1" means that this transition corresponds to an 
individual activity. The locations of "—1" and "1" indicate which local derivatives can enable 
the activity and what the derived local derivatives are, i.e. the pre and post local derivatives. In 
addition, the numbers of "— l"s and "l"s in each column are the same, because any transition in 
any component type corresponds to a unique pair of pre and post local derivatives. In fact, all 
this information is also stored in the labels of the activities. Therefore, with the transition rate 
functions defined in the next section to capture the timing information, a given PEPA model 
can be recovered from its activity matrix. 

Moreover, the pre and post activity matrix indicate the local derivatives which can fire a 
labelled activity and the derived local derivative after firing a labelled activity respectively. The 
modified activity matrix equals the difference between the pre and post activity matrices, i.e. 
C = C^°^* — C^'"''. Hereafter the terminology of activity matrix refers to the one in Definition|4j 
This definition embodies the transition or operation rule of a given PEPA model, with the 
exception of timing information. For a given PEPA model, each transition of the system results 
from the firing of an activity. Each optional result after enabling this activity corresponds to a 
relevant labelled activity, that is, corresponds to a column of the activity matrix. Conversely, 
each column of the activity matrix corresponding to a labelled activity, represents an activity and 
the chosen derived result after this activity is fired. So each column corresponds to a system 
transition. Therefore, we have the following proposition, which specifies the correspondence 
between system transitions and the columns of the activity matrix. 

Proposition 2. Each column of the activity matrix corresponds to a system transition and each 
transition can be represented by a column of the activity matrix. 

5 Transition Rate Function 

The structural information of any general PEPA model is captured in the activity matrix, which 
is constituted by all transition vectors. However, the duration of each transition has not yet 
been specified. This section defines transition rate functions for transition vectors or labelled 
activities to capture the timing information of PEPA models. 




+ 1 if Pi G pre(aj) 
otherwise. 



( +1 if Pi £ post(Q;j) 
I otherwise. 
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5.1 Model [2] continued 

Let us start from Model [2] again. As Table [s] shows, activity 7 in Model [2] is labelled as ^^a^^i 
and ^Qa-S'Qi^ Pqj. r^Pa^Pi ^ there are x[P3] instances of the component type P in the local 
derivative P3 in state x, each enabling the individual activity concurrently with the rate r^. So 
the rate of ^^s-s-A state x is /(x, = r-yx[P3]. Similarly, the rate for jQ^^Qi state 

X is r^'x[(52]- This is consistent with the definition of apparent rate in PEPA, which states that 
if there are replicated instances of a component enabling a transition (/, r), the apparent rate 
of the activity will be r x N. 

In Model [2] activity /3 is labelled as (3^'^^^'^ and /J^^-s-fs ^ ^o respectively reflect the corre- 
sponding two choices. According to the model definition, there is a flux of r^x(P2) into Pi from 
P2 after firing /3 in state x. So the transition rate function is defined as /(x, /3^2->A^ = r^x[P2]- 
Similarly, we can define /(x, /J^^^^b^ = r^x[P2]- These rate functions can be defined or in- 
terpreted in an alternative way. In state x, there are x[P2] instances that can fire /?. So the 
apparent rate of /? is (r^ -|- r|j)x[P2]- By the semantics of PEPA, the probabilities of choosing 

r' 

the outputs are ''f , and — r-r respectively. So the rate of the transition R^^^^^ is 

/(x,/3^^-^^) = ^^(r^ +r;,)x[P2] = r^^m (5-2) 

while the rate of the transition (3^'^^^^ is 

/(x,/3^^-^3) = +^;^)x[P2] = r-;,x[P2]. (5.3) 

In Model [2| a is a shared activity with three local rates: Tq,, and r^. The apparent rate of a 
in Pi is (r^-|-r^)x[Pi], while in Qi it is rQx[(5i]. According to the PEPA semantics, the apparent 
rate of a synchronised activity is the minimum of the apparent rates of the cooperating compo- 
nents. So the apparent rate of a as a synchronisation activity is min{(r^ -|- r^)x[Pi], rQ,x[(5i]}. 
After firing a. Pi becomes either P2 or P3, with the probabilities , ° „ and , " „ respectively. 
Simultaneously, Qi becomes Q2 with the probability 1. So the rate function of transition 
(Pi P2,Qi Q2), represented by /(x, a(^i-^-P2,Qi^Q2))^ is 



/(x,a 



r' 



(Pl^P2,Ql->Q2)) 

(5.4) 

-min{(r;, + r^)x[Pi],rax[gi]}. 



Similarly, 



/(x,a(^^^^3'«^^«^)) 

r" , „ (5.5) 

" — min{(r^ -Fr^)x[Pi],rQx[Qi]}. 



r'+r'' 



The above discussion about the simple example should help the reader to understand the 
definition of transition rate function for general PEPA models, which is presented in the next 
subsection. 



5.2 Definitions of transition rate function 

In a PEPA model, as we have mentioned, we may rewrite any U — — >■ V a.sU — > V, where r 
is denoted by r^^^ . The transition rate functions of general PEPA models are defined below. 
We first give the definition of the apparent rate of an activity in a local derivative. 
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Definition 5. (Apparent Rate of I in U) Suppose I is an activity of a PEPA model and U is 
a local derivative enabling I (i.e. U £ pre(/) ). Let post(?7, /) be the set of all the local derivatives 

derived from U by firing I, i.e. post(C/, /) = {V \ U — > V}. Let 

n{u)= '■y^''- (5-6) 

yGpost(i/,o 

The apparent rate of I in U in state x, denoted by n(x, U), is defined as 

ri{x,U)=x[U]ri{U). (5.7) 

The above definition is used to define the following transition rate function. 

Definition 6. (Transition Rate Function) Suppose I is an activity of a PEPA model and x 
denotes a state vector. 

1. If I is individual, then for each U ' — > V, the transition rate function of labelled activity 
[U^v state X is defined as 

/(x,F-^)=x[C/]r[^-^. (5.8) 

2. If I is synchronised, with pre(/) = {C/i,C/2r"" -.Uk}, then for each {Vi.,V2,--- ,Vk) in 
post(C/i, /) X post(C/2, X • • • X post([/fc, I), let w = (Ui — )• Vi, U2 ^ V2, ■ ■ ■ ,Uk ^ V^). 
Then the transition rate function of labelled activity l^ in state x is defined as 

fi^,n=[U\n] min {rKx,C/0}, 
where n(x, Ui) = x[[/j]r/([/j) is the apparent rate of I in Ui in state x. So 

/(x,r)= fn^Vv) f^^jmhm}. (5.9) 



\i=l 



Remark 1. Definition^ accommodates the passive or unspecified rate T. If some rf^^ are T, 
then the relevant calculation in the rate functions (5.8) and (5.9) can be made according to the 
inequalities and equations that define the comparison and manipulation of unspecified activity 
rates (see J^). Moreover, we assume that • T = 0. So the terms such as "minj AT, ri?} " are 
interpreted as |7iy ." 

. r ,^ ( rB, A> 0, 

mm{AT,r5} = | ^ ' ^ ^ 

The definition of the transition rate function is consistent with the semantics of PEPA: 

Proposition 3. The transition rate function in Definition^ is consistent with the operational 
semantics of PEPA. 

The proof is easy and omitted here. Since both the structural and timing information has 
been captured in the defined numerical representation schema, PEPA models can be therefore 
recovered from its representation schema. In addition, it is also easy to find that the transition 
rate function has the following homogenous property. 

Proposition 4. The transition rate function /(x, /) satisfies that for any H > 0, IIf{x/H, I) = 
/(x,/). 

This property will identify the CTMCs underlying a PEPA model to be density dependent 
(see Theorem [4] in the next section). 
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5.3 Algorithm for deriving activity matrix and transition rate functions 

This subsection presents an algorithm for automatically deriving the activity matrix and tran- 
sition rate functions from any PEPA model, see Algorithm [T| The lines 3-12 of Algorithm [T] 
deal with individual activities while lines 13 — 32 deal with shared activities. The calculation 
methods in this algorithm are essentially the embodiment of the definitions of labelled activity 
and apparent rate as well as transition rate function. So we do not give more introduction to 
this algorithm. 

6 Computational approaches for PEPA 

As a model being represented numerically, efficient techniques such as stochastic simulation and 
fluid approximation can be can be directly utilised to analyse the model. This section briefly 
introduces these approaches as well as technical foundations for employing them in the context 
of PEPA. 

6.1 Place/Transition structure in PEPA models 

Whilst the focus of stochastic process algebras has understandably been primarily quantitative 
analysis, qualitative analysis can also provide valuable insight into the behaviour of a system. 
In contrast, in Petri net modelling there are well-established techniques of structural analysis 
|12| - [T4] . This subsection shows how the new representation schema helps to manifest the P/T 
structure underlying PEPA models, and makes it possible to readily adapt structural analysis 
techniques for Petri nets to PEPA. First, the relevant definitions are given below. 

Definition 7. (P/T net, Marking, P/T system, Q/j 

1. A Place/Transition net (P/T net) is a structure M = (P, T, Pre, Post) where: P and T 
are the sets of places and transitions respectively; Pre and Post are the \P\ x |r| sized, 
natural valued, incidence matrices. 

2. A marking is a vector m : P — )• N that assigns to each place of a P/T net a nonnegative 
integer (number of tokens). 

3. A F /T system is a pair S = {M, mo).' a net Af with an initial marking mo. 

By Definition [TJ it is easy and direct to verify 

Theorem 1. There is a P/T system underlying any PEPA model, that is (A/", mo), where mo 
is the starting state; J\f = (T>,Aig^i)gi,C^^^,C^^^^^ is P/T net: where T> is the set of local 

derivatives, ^^^^g/ is the labelled activity set; qP^^ and qP^^^ are the pre and post activity 
matrices respectively. Moreover, each state m of the PEPA model is a marking. 

Based on the P/T structure underlying PEPA models and the theories developed for P/T 
nets, several powerful techniques and approaches for structural analysis of PEPA were estab- 
lished in [5] . For instance, the authors gave a method of deriving and storing the state space 
which avoids the problems associated with populations of components, and an approach to find 
invariants which can be used to qualitatively reason about systems. Moreover, a structure-based 
deadlock- checking algorithm was proposed, which can avoid the state space explosion problem. 
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Algorithm 1 Derive activity matrix and transition rate functions from a general PEPA model 
•^labei = 0; f is the set of all local derivatives 
for all activity / € ^ do 

if / is an independent activity then 
for all local derivatives U,V e V do 

ifU^V then 

^labei = ^labei U {l^^^} // Label I as 1^^^ 

1 1 Form a corresponding column of the activity matrix and the rate function 

r -1, d = c/ 

0, otherwise 
/(x,Z^^^) = rx[C/] 
end if 
end for 
end if 

if Hs a synchronised activity then 
pre(0 = 0, post(C/, l) = $,\/U eV 
for all local derivatives U,V e V do 

ifU^V then 

pre(0 = pre(Z) U {U} 
post(C/, I) = post(C/, I) U {V} 

end if 

end for 

Denote pre(/) = {pre(/)[l], pre(Z)[2], • • • , pre(Z)[A;]}, where k = #pre(/) 
for i = 1 . . .k do 

n(pre(OH) = E rr^')W-^ 

Vepost{pre{l)[i],l) 

end for 

K{1) = post(pre(Z)[l],0 x post(pre(/)[2], Z) x • • • x post(pre(/)[A;], Z) 
for all {Vi,V2,--- ,Vk) G K{1) do 

w = (pre(/)[l] ^ Vi,pre(0[2] ^ ^2, • • • ,pre(/)[A;] ^ Vk) 
^labei = ^labei U {P} // Label I as r 

/ / Form a column of Ma and the rate function corresponding to l'^ 

-1, d G pre(^) 
Ma{d,n = { 1, de{Vi,V2,--- ,Vk} 

0, otherwise 

/ k r,re(m^Vi \ 

= (n4^^(oM) j .„r„{'-'(p«(')W)xipre(owi} 

31: end for 
32: end if 
33: end for 

34: Output ^labei; M„; /(x, Z) (VZ G ^label)- 
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6.2 Stochastic simulation of PEPA models 



By solving the global balance equations associated with the infinitesimal generator of the CTMC 
underlying a PEPA model, the steady-state probability distribution can be obtained, from which 
performance measures can be derived. According to the original definition of the PEPA language 
in which each instance of the same component type is considered distinctly, the size of the 
state space of this original CTMC may increase exponentially with the number of components. 
By adopting the numerical vector form to represent the system state which results in the 
aggregated CTMC, the size of the state space can thus be significantly reduced, as Proposition [l] 
shows, together with the computational complexity of deriving the performance by solving the 
corresponding global balance equations since, the dimension of the infinitesimal generator matrix 
is the square of the size of the state space. 

Unless otherwise stated, hereafter the CTMC underlying a PEPA model refers to the aggre- 
gated CTMC, and the state of a model or a system is considered in the sense of aggregation. If 
the size of the state space is too large, it is not feasible to calculate the steady-state distribution 
and thus to get a performance measure i?, which is usually expressed as R = J2seS P(^)'^i^) 
where p and vr defined on the state space S are the reward function and the steady-state prob- 
ability distribution respectively. An alternative widely-used way to obtained performance is 
stochastic simulation. 

As discussed previously, a transition between states, namely from x to x + / , is represented 
by a transition vector / corresponding to the labelled activity I (for convenience, hereafter each 
pair of transition vectors and corresponding labelled activities shares the same notation). The 
rate of the transition / in state x is specified by the transition rate function /(x, /). That is, 

(Z,/(x,i)) , ^. . . . , . ,. r, . 

X — > x + 1. Given a starting state xq, the transition chain corresponding to a hrmg sequence 
^0, ^1, • • • ) ^ • • • is 

Xo — > xo + /o — > (xo + /o) + ti — > 

• • • > X 7- X -|- / > ■ ■ ■ . 

The above sequence can be considered as one path or realisation of a simulation of the aggregated 
CTMC, if the enabled activity at each state is chosen stochastically, i.e. is chosen through the 
approach of sampling. After a long time, the steady-state of the system is assumed to be 
achieved. Hence the average performance R = X^seS /^(^)^(^) be calculated. 

As one benefit of our numerical representation schema, it provides a good platform for 
directly and conveniently simulating the CTMC for PEPA, see Algorithm [2j In Algorithm [2| 
the states of a PEPA model are represented as numerical vector forms, and the rates between 
those states are specified by the transition rate functions which only depend on the transition 
type (i.e. labelled activity) and the current state. In this algorithm, the generated time r in each 

jution with the mean 

is in fact expressing: 



iteration can be regarded as having been drawn from an exponential distri 

, page 38). That is. Line 9 in Algorithm 



/(x) 



(see Example 2.3 in 
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"generate r from an exponential distribution with the mean ". Line 10 determines which 

/(x) 

transition will be chosen, and consequently determines the next state that the system will 
transition into. Therefore, this algorithm is essentially to simulate the CTMC underlying a 
PEPA model. 

The choices for stopping the algorithm include a given large time, or the absolute or relative 
error of two continued iterations being small enough (since the output performance converges 



^In practise, in order to decrease the computational cost, we should not calculate the performance until after 
a warm up period so that the effects of the initial state bias can be considered to be negligible. 
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Algorithm 2 Simulation algoritlim for deriving general performance measures from PEPA 

model 

1: / /Initialisation 

2: starting state x; labelled activity set ^labei = {hjhr" Jm}] activity matrix; transition 

rate function / 
3: reward function p; PerMeasure = 
4: while stop condition not satisfied do 
5: //Sampling 

6: compute the transition rate function f{x,lj), j = 1, 2, • • • , m 

7: /(x) = Er=l/(^'^.) 

generate uniform random numbers ri,r2 on [0, 1] 

compute r = — (l//(x)) Inri 

find fi such that Ylj=i /(x, Ij) < ?'2/(x) < Ylj=i /(x, Ij) 
I /Updating 

PerMeasure = PerMeasure + p(x) x r // Accumulate performance measur^ 
t = t + r // Accumulate time 
X = X + // Update state vector of system 
end while 

Output performance: P^-'Measure 
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as time goes to infinity as the following Theorem [3] states). Now we prove the convergence of 
the performance calculated using the algorithm. We need the following theorem. 

Theorem 2. (Theorem 3.8.1, fW) If X{t) is an irreducible and positive recurrent CTMC with 
the state space S and the unique invariant distribution n, then 



Pr ^ l{x,=s}dz 7r(s) as t ^ oo^ =1. (6.10) 



Moreover, for any bounded function p : S" — t- M, we have 



Pr ^ p{X^)dz E[p{X)] ast-^oo^ =1. (6.11) 

where E[p{X)] = E^^gPi^M^)- 

Here is our conclusion (we assume the CTMCs underlying PEPA models to be irreducible 
and positive recurrent): 

Theorem 3. The performance measure calculated according to Algorithm^ converges as time 

goes to infinity, that is, limt^^ PerMeasure ^ E[p{X)]. 

Proof. Assume that n — 1 iterations have been finished and the time has accumulated to tn-i- 
Suppose the current one is the n-th iteration and r is the generated time in this iteration. After 
the n-th iteration is finished, the accumulated time will be updated to t„ = + r. During 
the r time interval, the simulated CTMC stays in the state x, that is, X^ = x, z £ [t„_i,t„). 
So, 

p{x) X T 



p(x) / dz= p{x)dz = / p{X,)dz. 

^ i-n — l i-n — l ^n — 1 
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Therefore, after this n-th iteration, PerMeasure will be accumulated to = J"q" p{Xz)dz and 

PerMeasure 1 /"*" 

tn tri Jo 



p{X,)dz. 



According to Theorem 



1 

— / p{Xz)dz tends to E[p{X)] as tn tends to infinity. So the 
j_, Jq 

performance obtained through Algorithm ^ converges to E[p{X)] as the simulation time goes 
to infinity. □ 



Performance metrics, such as activity throughput of an activity and capacity utilisation of 
a local derivative that are discussed in j6], can be derived through this algorithm by choosing 
appropriate reward functions. 

6.3 Fluid approximation of PEPA models 

The weakness of the simulation method is its high computational cost, which makes it not 
suitable for real-time performance monitoring or prediction. Recently, a novel approach to get 
performance measures from PEPA models has been proposed in [4j and subsequently expanded 
in [T] and [17| , making a continuous state space approximation as a set of ordinary differential 
equations (ODEs). In this section, we present a mapping semantics for this approaches, which 
is based on the numerical representation schema. In addition, a theoretical justification of this 
approach, mainly in terms of its consistency with the CTMCs, will be discussed. 

In our representation schema, the transition rate /(x, /) refiects the intensity of the transition 
from state x to state x + /. The state space is inherently discrete with the entries within the 
numerical vector form always being non-negative integers and always being incremented or 
decremented in steps of one. As pointed out in [4], when the numbers of components are 
large these steps are relatively small and we can approximate the behaviour by considering the 
movement between states to be continuous, rather than occurring in discontinuous jumps. In 
fact, let us consider the evolution of the numerical state vector. Denote the state at time t by 
x(t). In a short time At, the change to the vector x(t) will be 

x(-,t + At)-x(-,t) = At lf(.^{;t),l). 

Dividing by At and taking the limit. At — )• 0, we obtain a set of ODEs: 

^= Yl ^/(^'O- (6.12) 

Once the activity matrix and the transition rate functions are generated, the ODEs are 
immediately available. All of them can be obtained automatically by Algorithm [TJ 

For an arbitrary CTMC, the evolution of probabilities distributed on each state can be 
described using linear ODEs (see [Is], page 52). For example, for the aggregated CTMC un- 
derlying a PEPA model, the corresponding differential equations are 

^ = Q^vr, (6.13) 

where each entry of 7r(t) represents the probability of the system being in each state at time t, 
and Q is an infinitesimal generator matrix corresponding to the CTMC. Clearly, the dimension 
of the coefficient matrix Q is the square of the size of the state space, which increases with the 
number of components. 
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The scale of (6.13), i.e. the number of the ODEs, depends on the size of the state space, 



so it suffers from the state sapce explosion problem. In contrast, the ODEs (6.12) reflect the 



evolution of the population of the components in each local derivative, so the scale is only 
determined by the number of local derivatives and is unaffected by the size of the state space. 
Therefore, it avoids the explosion problem. But the price paid is that the ODEs (6.12) are 



generally nonlinear due to synchronisations, whereas (6.13) is linear. 



This paper emphasises the consistency between the fluid approximation and the aggregated 
CTMC. Obviously, the CTMC depends on the starting state of the given PEPA model. By 
altering the population of components presented in the model, which can be done by varying 
the initial states, we may get a sequence of aggregated CTMCs. Moreover, the homogenous 
property that the transition rate function satisfies, indicated in Proposition [4j identifies the 
aggregated CTMC to be density dependent. 

Definition 8. ( [191). A family of CTMCs {X„}„ is called density dependent if and only if 
there exists a continuous function /(x, /), x G M'^, / e Z'^, such that the infinitesimal generators 
of Xn are given by: 



Q 



(n) 



n/(x/n,0, ly^O, 



where 9xx+/ denotes an entry of the infinitesimal generator of Xn, x a numerical state vector 
and I a transition vector. 

This allows us to immediately conclude the following conclusion: 

Theorem 4. ( f^\11^) Let {Xn} he a sequence of aggregated CTMCs generated from a given 
PEPA model (by scaling the initial state), then {Xn} is density dependent. 

Since both ODEs and density dependent CTMCs can be derived from the same PEPA model 
through the same activity matrix and transition rate functions, it is natural to believe some kind 



of consistency between them. In fact, according to Kurtz's theorem 20 , the complete solution 



of some ODEs can be the limit of a sequence of Markov chains. Such consistency in the context 
of PEPA has been previously illustrated for a particular PEPA model in ^21j , and subsequently 
generalised to general models in [7] and 17 . The result presented below is extracted from 17 



in which the convergence is in the sense of almost surely rather than probabilistically as in 
and [7]. 
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Theorem 5. ( flTV) Let x(t) be the solution of the ODEs (6.12) derived from a given PEPA 



model with initial condition xq, and let {Xn{t)} be the density dependent CTMCs with X„(0) 
nxQ. Then for any t > 0, 



lim sup — x(n)|| =0 a.s. 

u<t 



(6.14) 



This theorem justifies the fluid approximation by manifesting the consistency between this 
approach and the corresponding CTMCs for a general PEPA model. Furthermore, if there is 



no synchronisation contained in the model then the derived ODEs (6.12) becomes linear, and 



(6.13) and (6.12) coincide except for a constant factor. Moreover, the fundamental results on 



the fluid approximation of PEPA models such as the existence, uniqueness, boundedness and 
nonnegativeness of the ODEs' solution, as well as the solution's asymptotic behaviour, have 
been obtained. In particular, the convergence of the ODEs' solution as time tends to infinity, 
has been proved under a condition, which is revealed to relate to some famous constants of 
Markov chains such as the spectral gap and the Log-Sobolev constant. For more details about 
these stories, please refer to |22j . As for performance derivation via this approach, please see 
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7 Conclusions 



In this paper we have demonstrated a schema, which bridges the syntactic and numerical repre- 
sentation, as well as the local definition and global analysis for a PEPA model. Computational 
approaches and associated algorithms developed based on the schema have been presented, 
which can help to relieve the state space explosion problem for large scale models. For other 
stochastic process algebras, similar numerical representation schema can be established and 
expected to benefit relevant performance modelling. 
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